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Abstract. We study a self-similar circulation model for 
protostellar bipolar outflows. The model is axisymmct- 
ric and stationary, and now includes Poynting flux. Com- 
pared to an earlier version of the model, this addition 
produces faster and more collimated outflows. Moreover 
the luminosity needed for the radiative heating is smaller. 
The solutions are developed within the context of r-self- 
similarity, which is a separated type of solution wherein a 
power of r multiplies an unknown function of 0. For out- 
flows surrounding a fixed point mass the velocity, density 
and magnetic field respectively scale with spherical radius 
rasoa r -1 / 2 , p ex r 2 " -1 / 2 and B ex r Q ~ 3 / 4 . The param- 
eter a must be larger than —1/2 and smaller than or equal 
to 1/4. We obtain the ^-dependence of all flow quantities. 
The solutions are characterized by a together with the 
scaled temperature parameter O that imposes the neces- 
sary heat transfer. The model has been applied to both 
low and high mass protostars. Monte Carlo methods have 
been used to explore systematically the parameter space. 
An inflow/outflow pattern including collimation of high 
speed material and an infalling toroidal disc arises natu- 
rally. The disc shape depends on the imposed heating, but 
it is naturally Keplerian given the central point mass. Out- 
flows can have large opening angles, that increase when 
magnetic field weakens. Massive protostars produce faster 
but less collimated outflows than less massive protostars. 
The model is now at a stage where synthetic CO spectra 
reproduce very well the observational features. The results 
strengthen the idea that radiative heating and Poynting 
flux are ultimately the energy sources driving the outflow. 
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1. Introduction 

General outflows and collimated jets are thought to be 



cent observations (Donati et al. 1997, Guenther & Emer- 
son 1996[ ) suggest the relevance of magnetic fields to star 
formation, first suggested by Mestel & Spitzer (1956). 
Tomisaka ( |1998| ) has studied numerically the dynamical 
collapse of magnetized molecular cloud cores from the 
runaway cloud collapse phase to the central point mass 
accretion phase. He finds that the evolution of the cloud 
contracting under its self-gravity is well expressed by a 
self-similar solution. Moreover inflow-outflow circulation 
appears as a natural consequence of the initial configura- 
tion. His results support the magnetized self-similar mod- 
els as recently presented by Fiege and Henriks en (1996 a,b, 
hereafter FH1, FH2) following Henrik sen ( |1993| , |l994j ), 
and Henriksen and Valls-Gabaud (1994, hereafter HV). In 
those models, the self-similar gravitationally driven con- 
vective circulation in a heated quadrupolar protostar en- 
velope is solved rigorously. Moreover recent observations 
(Greaves & Holland 1998) of the polarized dust emission 
from six star-forming cloud cores have revealed, for the 
first time, some very large twists in the magnetic field. 
This is consistent with FH1 and HV where circumstel- 
lar magnetic fields follow large-scale streamlines. This re- 
mains very nearly the case in the present model. 

The self-similar models regard the molecular outflow 
as a natural consequence of the circulation established by 
the collapse of the pre-stellar cloud. The models describe 
an outflow velocity that increases toward the axis of rota- 
tion, a convective pattern of inf all/ outflow, self- consistent 
axial collimating magnetic fields and rotation, and "cored 
apple" type distributions of circum-proto stellar gas (Andre 
et al. 1199$ . 

Such models may imply a natural connection between 
the fast ionized jets seen near the polar axis of the wind, 
and the slower and less-collimated molecular outflows that 
surround them. However velocities obtained by FH1 and 
HV are smaller than some observed velocities (e.g. EHV 
outflow (V v 100 km s" 1 ), Backfiller et al. |1990| ) unless 
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nant and the assumptions of this type of model are liable 
to break down. Moreover the luminosity needed for ra- 
diative heating in FH1 is too large in order to drive high- 
velocity outflows if only dust opacity is used. In the present 
work we show that the Poynting flux increases both the ve- 
locity and collimation of the modeled outflows by helping 
to transport mass and energy from the eguatorial regions. 
This is much as has been argued for some time by other au- 



thors (e.g. Pudritz and Norman 1986, Ouyed and Pudritz 



1997 ), but our models are globally consistent in space. 



Self-similar models can not of course be globally consis- 
tent in time even if they are not stationary (we have found 
such models and will report them elsewhere), since they 
are ignorant of initial conditions. And in fact they must 
also be 'intermediate' in space since there are inevitably 
small regions near the equator, the axis and at the cen- 
tre which are excluded from the domain of self-similarity. 
Thus, we derive globally self-consistent models of bipolar 
outflows and infall accretion, but specifically exclude the 
regions dominated by the disc, jet, and protostar. 

Qualitatively we may understand the quadrupolar 
global circulation in the following fashion. As material falls 
towards the central object, it is gradually slowed down 
by increasing radial pressure gradients due to the rising 
temperature, density, and magnetic field encountered by 
material near the central object. This pressure barrier, 
with the help of centrifugal forces, deflects and acceler- 
ates much of the infalling matter into an outflow along 
the axis of symmetry. The outflow velocity is naturally the 
highest near the axis of rotation since the pressure gradi- 
ents are strongest there. Mathematically the quadrupolar 
model can be seen as a form of instability of the spheri- 
cal Bondi accretion flow, wherein rotation, magnetic fields 
and anisotropic heating transform a central nodal singu- 
larity into a saddle point singularity. 

The pressure gradients also act to collimate the flow, 
forming consistently the now traditional de Laval nozzle 



(e.g. Blandford and Rees |l97j, Konigl |l982[ ). We find how- 
ever that the most convincing flows are everywhere super- 
Alfvenic, in agreement with FH1, but the nozzle may still 
be critical in terms of the fast magnetosonic speed as dis- 
cussed in FH1 and below. We note that one can expect 
conical shocks near the rotational axis, which would seem 
to be necessary to explain the shocked molecular hydrogen 
emission often observed. 

The paper is organized as follows: in §2 we introduce 
the model and its approximations, and discuss their conse- 
quences; numerical solutions are presented in §3 for virial- 
isothermal and radiative models describing outflows from 
low and high mass protostars. In §4 we study the ensuing 
synthetic radio maps; this section is followed by discus- 
sions of the results in §5, and conclusions are presented in 



2. The Analytical Model 

The main difference between hydrodynamic and MHD 
outflows is the low asymptotic Mach number in the MHD 
case. For protostellar outflows, the magnetosonic speed is 
about 100 km s~ l at 10 5 AU and therefore far exceeds 
the sonic speed. It is not a priori obvious that the internal 
Alfvenic Mach number of the jet or outflow is small. How- 
ever most of our models are super- Alfvenic everywhere. 

The plasma is supposed to be perfectly conducting and 
therefore the flow is governed by the basic equations of 
ideal MHD without taking into account resistivity or am- 
bipolar diffusion. In order to make the system tractable 
we assume axisymmetric flow so that d/dcj) — and all 
flow variables are functions only of r and 9. 



Kudoh and Shibata ( 1997b ) have performed time de- 
pendent one-dimensional (1.5-dimensional) magnetohy- 
drodynamic numerical simulations of astrophysical jets 
that are magnetically driven from Keplerian discs. They 
have found that the jets, which are ejected from the disc, 
have the same properties as the steady magneticall y drive n 
jets they had found before (Kudoh & Shibata 1997a). 
Their numerical results suggest that a steady model, as 
we are assuming in the present work (i.e. d/dt = 0), is a 
good first approximation to the time-dependent problem 
on time-scales short compared to that required to dissi- 
pate the circumstellar material. Time independence will 
be relaxed in future work. 

Accretion discs near protostars are probably heavily 
convective and therefore prone to dynamo action (Bran- 
denburg et al. 1995, FH1). Since the ensuing magnetic 
fields are loaded with disc plasma, currents are allowed 
to circulate between the disc surface and the wind region 
above the disc. For rapidly rotating discs, this dynamical 
system may evolve naturally into a quadrupolar structure 
(Camenzind 1990J , Khanna and Camenzind 1996). The 
strong differential rotation of accretion discs is responsible 
for the excitation of quadrupolar modes in discs. Conse- 
quently, in the present model, magnetic field and stream- 
lines are required to be quadrupolar in the poloidal plane. 
A moderate amount of heating supplied near the turn- 
ing point allows the outflow to possess finite velocities at 
infinity. 

These considerations have led us to use the set of equa- 
tions of steady, axisymmetric, ideal MHD and to seek so- 
lutions with a quadrupolar geometry. Even though the 
poloidal components of the magnetic and velocity fields 
have to be parallel, the toroidal components need not 
share the same constant of proportionality (Henriksen 
1996). This permits a poloidal conservative electric field 
to exist in the inertial frame, and so admits steady Poynt- 
ing flux driving. One needs then to introduce an electric 
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transfer (to zero order in v/c) and a Kramer's-type law 
for the Rosseland opacity are added as in FH1. We also 
neglect the self-gravity of the protostellar material by as- 
suming that the gravitational field is dominated by the 
central mass that is assumed to be an external fixed pa- 
rameter in our model. Henceforth the model only applies 
to a state where the protostellar core is already formed. 

It is unlikely that this model can be extended to in- 
clude the optical jets since these must originate very near 
to the central object where the model may not apply in its 
steady form. Nevertheless, the inclusion of a jet model i n 
the axial region, such as given by lery et al.( 1998 , 1999b|) , 
could remedy this lack in order to make a more global 
model for protostellar objects. 

2.1. Scales and Self-similar Laws 

Our model is developed within the context of powe r law 
self-similar models as developed by Henriksen ( 1993 ), and 
in HV and FH1. A model with such symmetry was al- 
ready used by Bardeen & Berger ( 1978J ) for galactic winds, 
although the present development has been independent. 
Parallel developments have been largely occupied with the 
study of winds from an establi shed accretion disc (Bland- 
ford & Payne |l982|, Konigl |l989j, Ferreira & Pelletier 



1993| , Rosso fc Pe lletier |1996| Ferreir a |1997| , Contopoulos 
fc Lo velace 1994 , Con topou los 1995 , Sauty & Tsinganos 
1994 Tsinganos et al. |l996| , Vlahakis & Tsinganos |l99^ ) . 



Our use of this model in order to study inflow and outflow 
as part of the global circulation dynamics around proto- 
stars seems unique. The form used corresponds to an ex- 
ample of 'incomplete self-similarity' in the classification 



of Barenblatt (1996), but fits into the general scheme of 
self-similarity advocated by Carter and Henriksen (1991) 
with r as the direction of the Lie self-similar motion. We 
work in spherical polar coordinates r, 9, <f> centered on the 
point mass M and having the polar axis directed along the 
mean angular momentum of the surrounding material. 

The self-similar symmetry is identical to that used in 
FH1 except that two quantities y^ and y p defined such 
that i/AKJIyp = vpjiBpj \J^-Kp) (where (3 indicates the 
appropriate component) replace y in FH1. The power laws 
of the self-similar symmetry are determined, up to a single 
parameter a, if we assume that the local gravitational field 
is dominated by a fixed central mass. In terms of a fiducial 
radial distance, r , the self-similar symmetry is sought as a 
function of two scale invariants, r/r and 9, in a separated 
power-law form. The equations of radiative MHD require 
the following radial scaling relations for the variables: 
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In these equations the microscopic constants are repre- 
sented by k for Boltzmann's constant, p, for the mean 
atomic weight, run for the mass of the hydrogen atom. 
The self-similar index a is imposed as a parameter of the 
solution, but as in FH1 it is constrained to lie in the range 
1/4 > a > —1/2. In the last equation, the index a/ is a 
measure of the loss (if negative) or gain in radiation energy 
as a function of radial distance. 

Whenever we take account of radiative diffusion ex- 
plicitly, as well as for the definition of the luminosity, we 
will use the same formulation as in FH1. The simplest 
approach however is to assume that the temperature is 
some fraction of the 'virialized' value. That is, according 
to relation (g), that O takes some constant value < 1. By 
using the first law and assuming an ideal gas (so that with 
the above relations for p and p we have P{9) = jjlQ ) this 
assumption implies that on a sphere the specific entropy 
s varies according to 

kT 
Tds = -- dlnu. (8) 

fim H 

Consequently we are explicitly adding or subtracting heat 
from the system on a spherical surface according to the 
sign of — dln/x. In our models this normally increases to- 
wards the rotational axis. Thus, at least implicitly, we al- 
ways employ some form of heat transfer, presumably ra- 
diative, to the gas over the poles relative to the equatorial 
material. 

The same procedure for variations in the radial direc- 
tion yields 
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2 dr 
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Energy is therefore lost to the material with increasing 
radius provided that the quantity in the square bracket 



is positive. This requires a — 1/4 > 



which is 



' 2(7-1) : 

true for reasonable ratios of specific heats throughout our 
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of the dependence on F raa ; this is not the same as requir- 
ing an optical depth of order unity. Nevertheless in order 
for the radiative transfer by diffusion to be plausible, this 
latter condition must also be satisfied. 

In either case the issue depends on the nature of the 
opacity k. We follow FH1 in supposing that this is pre- 
dominantly a dust opacity (at least at some distance from 
the 'star') which is taken in Kramer's form 



Ko P a T b 



(10) 



A fit to a dust model in FH1 yields a — 0, and b ~ 2 and 
k = 6.9 x 10~ e cgs. We can calculate the mean optical 
depth between an inner radius r* and the fiducial radius 
r , namely < r >= j- J ° npdr dtt, using this opacity 
and the self-similar scaling relations. We find 

</i6 2 
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where we have defined the fiducial radius as 

1/4 
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G k 



M 3 / 4 , 



which is numerically 

r w U00(M/M Q ) 3/4 AU. 



(11) 



(12) 



(13) 



Thus we obtain a fiducial radius that is characteristic of 
the bipolar outflow sources and is a convenient unit for our 
discussions. A similar scale (to withi n a factor 5) was pre- 
viously deduced in Henriksen ( 1993 ) based on the source 
luminosity and fundamental constants (including the ra- 
diation constant). We note that the external opacity from 
oo to r is given simply by < t > ex t—< M@ 2 > /|2a— 3/2| 
and the 'photosphere' of the cloud is therefore found where 

Text ~ 2/3. 

The corresponding temperature at the fiducial radius 
is given by 



r« 164Q (M/MqY^K. 



(14) 



We observe that if we identify our radiative fiducial scale 
r with an empirical scale r em f rom < [i > /r 2 m — 
2 x 1(T 9±1 ' 



1998), we can infer a V a e = 



l (AU)- 2 (Richer et al. 
relation between the mean density parameter < fi > and 
the central mass in the form 



< n >- 



3.38 x 1CT 3±1 (M/M e ) 3/2 . 



(15) 



This is a useful way to relate the key model parameter 
jjl to the physical mass and luminosity. We note however 
that the optical depth is not likely to approach unity for 
reasonable temperatures, except for exceptionally massive 
objects. But this does not prevent radiative heating from 
playing an important dynamical role near the protostar. 
In fact using r , as above, and the definition of ra- 



The radiative force per unit mass is given by 

nF rad /c « 10- 13 (Af/M ) 1/8 (r/r ) (Q '" 4) fern s~ 2 , (17) 

which must be less than the Eddington limit of GM/r 2 = 
3.5 x 10~ 7 (M/M Q ) (r /r) cm s~ 2 , for consistency. These 
are helpful considerations for our radiative models below 
that always verify this condition a posteriori. 

2.2. The Equations 

We use the self-similar forms in the usual set of ideal MHD 
equations together with the radiative diffusion equation 
when applicable. The corresponding system of equations 
is reported in Appendix A, and the ensuing first integrals 
in Appendix B. As noted above, in the present model, 
there is no constraint requiring parallelism between veloc- 
ity and magnetic field as used in FH1. This allows us to 
discern clearly the effects of a non-zero Poynting flux. The 
self-similar variable directly related to magnetic field y has 
two components that are not equal in the present model. 
These quantities are either projected in the poloidal plane 
or correspond to the toroidal components and are respec- 
tively defined by 



VpA°) 
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(18) 



Consequently the system deals with two different compo- 
nents of the Alfvenic Mach number M ap and M at p defined 
by: 



M, 



ap.a,(p 
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(19) 



In this model, since r and 4> respectively correspond 
to the directions of self-similarity and axisymmctry, only 
waves that propagate along 9 in the poloidal plane can pre- 
serve these two symmetries. Therefore the relevant Alfven 
mode propagates in the direction 9 in the poloidal plane 
with a phase velocity: 



B e 



(20) 



Moreover we define an Alfvenic point where the total 
poloidal flow velocity is equal to this value in magnitude 
and direction (i.e. u r = 0), which by (ttq) is equivalent to: 



M ap = 1 = 4:TTfiy 2 . 



(21) 



One can then define a critical scaled density corresponding 
to the Alfvenic scaled density 
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8 with phase velocities that satisfy the quartic (Sauty and 
Tsinganos |1994| ) 



{v., f ) 4 g - (v s/f )l (v 2 + c 2 s ) + c 2 (V a ) 2 e = 0, 



(23) 



where C 2 is the isothermal sound speed. The condition 
that ug be equal to one of these wave speeds (where we 
expect a singular point in the family of solutions) can be 
expressed in terms of self-similar variables as: 



e 



2 2—2 
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and in physical units (cf FH1): 
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Equation (£3|) shows that whenever V a » C s the slow 
and fast speeds become V a $C s /V a and \/V 2 — V 2 g C 2 /V 2 
respectively. Thus a super- Alfvcnic flow as defined above 
can only encounter the fast speed and then only where 
V a > V a e- In a region where the sound speed domi- 
nates, the fast and slow speeds become \JC 2 — V 2 and 
V a e respectively. Either the sound speed or the Alfven 
speed is dominant according as 8 is larger or smaller than 
ul/M^ + ul/M^. 

As in Henriksen ( 1996| ) we can consider the energy 
balance in these models. A convenient form for the energy 
equation is 



V. P v( 
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where the Poynting flux vector is 
B.v 
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and $ = —GM/r. In self-similar variables the Poynting 
flux becomes 
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where the self-similar Poynting flux cr(9) can be defined 
in terms of the other variables as: 
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The energy equation becomes 
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2.3. The Zero Pressure Limit 

In the zero pressure limit one can use the self-similar en- 
ergy equation together with a r = (u r /ug)ag and u r /ug = 
d In r/d9 on a stream line to write a Bernoulli integral 

r 2a sin 6 [fiu e {u 2 /2 - 1) + a e ] = B, (31) 

where B is the Bernoulli constant on each stream line. In 



(32) 



view of Eq. 


B.2 we can write this as 
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where the second term evidently describes the Poynting 
driving (magnetic energy term) and the first term gives 
the kinetic and gravitational energies. Hence if B ^ 0, we 
require the expression contained in brackets on the RHS 
of this last equation to tend towards infinity both at the 
equator and on the polar axis since r(Q) — ► oo in these 
limits. This is an improbable outer boundary condition 
although physically realizable at some large but finite ra- 
dius (which might be a function of angle and so trace out 
an excluded region for the self-similar solution). It does 
permit the exchange between say a large driving Poynt- 
ing term near the equator and a large kinetic energy term 
near the axis, but the source of the 'free' magnetic energy 
is left undetermined. 

We turn to the special case B = 0, which requires this 
same factor to be zero everywhere on a streamline. In or- 
der that there be an exchange between the kinetic term 
and the Poynting term in the sum we must have either 
u 2 < 2 and y^ < y p or the reverse of both inequalities. 
In the first case we do not produce velocities greater than 
the escape velocity anywhere, while in the second case we 
require M ap < M a § everywhere on the streamline which 
also seems improbable. We are left then with the first pos- 
sibility, which includes the case where each term vanishes 
separately on the boundaries. This suggests that starting 
from a physically likely free-fall boundary condition, v is 
less than or equal to the escape velocity on every stream 
line in the absence of pressure, becoming equal again to 
the escape velocity at large distances. 

Thus without pressure effects one is unable to produce 
net outflows at infinity without appealing to a Poynting 
flux from an excluded equatorial region. However one does 
thus deflect and collimate the material (magnetically and 
centrifugally) in a lossless way, and so all energy added 
near the star may appear in the outflow at infinity rather 
than be lost to work against gravity. The fact that there is 
no net gain of kinetic energy from the magnetic field en- 
ergy on a stream line is peculiar to the case with B = 0. 
It is due to the fact that energy that is first gained by 
the magnetic field at the expense of gravity and rotation 
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3. Numerical Analysis 

3.1. The Numerical Method 

The numerical solutions in this work are obtained by initi- 
ating the integrations of the system of six equations from 
a chosen angle between the axis (9 = 0) and the equator 
(9 = §). Starting values are varied until the boundary con- 
ditions are met. Solutions can not strictly reach the axis 
and therefore are limited by a minimum angle Q m in- More- 
over ? bounds the solution only if we demand exact reflec- 
tion symmetry about the equator. Solutions are thus gen- 
erally defined on a wedge given by < 6 m i n < 9 < 9 max . 
In order to be able to start at will from either a sub- 
or super- Alfvenic point, a systematic search for singular 
points of the system has been conducted. This allows us to 
restrict a first guess for the typical values to use as input 
for the system. 

Several general criteria have been employed to fur- 
ther constrain the solutions. First the radial velocity u r 
must vanish once, to have inflow-outflow, and only once, 
to avoid a higher order convection pattern. We also re- 
quire ug and u<f, to vanish at both boundaries, since there 
must be no mass flux across them, and since only zero 
angular momentum material can reach the axis (which all 
material lying strictly in the equator does). In practice, 
the first of these two conditions is sufficient to impose the 
second one. The same constraints apply to the magnetic 
field components. We take ug to be negative in order to 
get a circulation pattern rising over the poles, and u^ to be 
positive as a convention. We require also that u r is either 
decreasing or constant with increasing angle and that the 
density and pressure attain their maxima close to 9 max . 

To find solutions satisfying all of the previous require- 
ments a Monte Carlo shooting method has been used 
starting from values corresponding to cither sub- or super- 
Alfvenic starting points. Solutions have been found using 
as starting points successively 9 m i n , 9 max and the stream- 
line turning point where the radial velocity vanishes. It 
has been found that the most convenient starting point 
is 9 m i n . Results of these investigations and representative 
cases are shown below. 



3.2. Low Mass Protostars 

Most of the present models for star formation deal only 
with low mass stars. It is then natural to begin with this 
type of object. The solution shown in Fig.H is represen- 
tative of this class in the virial-isothcrmal approximation. 
It shows the variations with angle of the six self-similar 
variables (related to velocity (u r , ug, u^,), density (/i and 
fi a ), magnetic field (y p , y^)) together with poloidal and 
toroidal Mach numbers (M ap , M a( f,) and poloidal compo- 
nent nf the Pmmtinc flnv (\*?_1 Tn this cnhitinn we have 



15 


-i ' 1 ' 1 ' r- 




' 1 ' 1 ' r- 


10 


-\ U 


-0.2 


u e/ 


5 




-0.4 


/ - 







-0.6 

0.0015 
0.0010 
0.0005 

10 4 

io 3 


^~~~~-------^__/ 




A u * 


_—-^ 


1 
0.5 




10* 
l(f 


y P - 


y, i 


10" 
10 1 




10" 

io 1 


V__ j 






10 3 
10 2 


— M - 

ap 
^V a"t» 


I A Sp 


io 1 


- ^-^_^ j 


io 4 


- v^ ^y v___^ 


10" 


^ -i _«rf< 




■ 



0.5 



1.5 



0.5 



1.5 



Fig. 1. Variations with angle of the self-similar variables 
(related to velocity (u r , ug, u^), density (/j, and fj, a ), mag- 
netic field (y p , ytj,) and Mach number (M op , M a $) for a 
virial-isothermal low mass solution with a = —0.3 and 
= 0.41. The poloidal component S p of the Poynting 
flux is plotted in the lower right panel. Mwl M Q . 



both vanish at the polar and equatorial boundaries, as do 
the magnetic field components (y is diverging). The max- 
ima for rotation velocity and Poynting flux, correspond- 
ing to the minimum in ug, occur where the radial velocity 
changes signs. This point defines the total opening angle 
of the outflow which in this model is of the order of 60°. 
The radial component of the Poynting flux S r changes 
signs as u r does, but its poloidal component defined by 
S p = y/S^+Sg remains positive by definition and shows 
only a slight dip near maximum. 

The high velocity outflow is of course much more nar- 
rowly collimated than the general outflow (u r > 0). The 
density fi is always superior to the Alfvenic critical density 
and both densities are almost equal where u r changes sign. 
Thus the flow is always super- Alfvenic as confirmed by the 
Mach number components in Fig.|l|. Since differences be- 
tween y p and y^ are not evident to the eye in this figure, 
we have plotted the poloidal component of Poynting flux 
as given by Eq. |29|. It shows the peak energy carried by 
the magnetic field to be at the turning point of the flow. 

Thus this calculation produces wide angle relatively 
slow outflows surrounding a fast component that is iden- 
tifiable with the EHV outflows or 'molecular jets'. We shall 
see below that inclusion of the Poynting flux allows more 
collimated and faster 'jets' surrounded by a slower wide 
angle outflow than is the case without it. Therefore the 
ma.enetic field nla.vs a, crucial role in the YSO outflow de- 
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Fig. 2. Streamlines of the flow in the poloidal plane for a — —0.3 and O = 0.41. They are integrated out to a radius 



of 4000 AU. Panel A shows a solution that almost fills full quadrant (from 9 
smaller opening angle and goes from 9 w to 1.33. 



Otowf). Solution in panel B has a 



size of about 10 4 AU (Adams et al. 1993, Terebey et al. o.4 



199?, Andre & Montmerle 1994). The size that we will use 



to present our solutions will be of this order. For example, 
two poloidal projections of the streamlines are represented 
in Fig.0 out to a radius of 4000 AU. On the left panel (case 
A), the previous solution is shown from (nearly) 9 = to 
(nearly) -|. In the panel on the right (case B) we show 
the result of varying the magnetic field and density from 
the previous values. The solution shown stops at 9 — 1.33 
and the empty region is outside of the domain of the so- 
lution. If the pressure and tangential velocity matching 
were done properly at the boundaries , a ro tating Bondi ac- 
cretion flow (Henriksen and Heaton 1975 ) might describe 
naturally a true accretion disc. 

The self-similar structure of the streamlines does not 
in fact strictly match an arbitrary boundary condition at 
infinity, but the calculated symmetry seems to be a nat- 
ural small scale limit of a more general self-gravitating 
circulation flow. Moreover the maximal meridional veloc- 
ity close to the axial region (see Fig.fil) shows the ten- 
dency for the flow to collimate cylindrically along the axis 
of rotation. This is in agreement with Lery et al 



1998 



1999b) where it is shown that cylindrical collimation is a 



general asymptotic behavior of magnetized outflows sur- 
rounded by an external pressure, due in the present case 
to the outer part of the molecular cloud or the interstellar 
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Fig. 3. View down the axis of streamlines in the axial 
region for the virial-isothermal solution (a = —0.3,0 = 
0.41). The streamline plotted with a solid line has the 
smallest angle with respect to the rotational axis. Other 
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Fig. 4. Virial-isothermal low mass solution: the angular 
dependence of magnetic field (in Gauss) at a distance of 
5 x 10 4 AU. Total magnetic field and r,0,(p components of 
B are represented, (a = -0.3,6 = 0.41,M « 1 M ) 



rial region and end in the axial region with small angles 
(9 = 0.005 to 0.15). In addition to the meridional behav- 
ior shown previously in Fig|2, each streamline makes a 
spiraling approach to the axis and then emerges in the 
form of an helix wrapped about the axis of symmetry. 
The closer in angle to the equator the streamline starts, 
the larger the angle of rotation it makes around the axis, 
the closer it gets to the central mass, and the nearer to 
the axis it asymptotically emerges. On the other hand, for 
streamlines beginning well above the equator, rotation is 
very nearly negligible and the path almost remains in a 
poloidal plane. 

Indirect measurements of the magnetic field in proto- 
stars are presently available (inferred from the radio flux 
from gyro-synchrotron emission by Ray et al. 1997, and 



Hughes 1999). Such measurements can provide a critical 
test of models. In Fig.[| the angular dependences of the 
components of the magnetic field are shown in Gauss at a 
distance of 5 x 10 4 AU for the same illustrative example of 
a low mass protostar. Magnetic fields of the order of a mil- 
ligauss or less are obtained with a peak around 1 mGauss 
that coincides in angle with the radial turning point of the 
flow. As required by the quadrupolar geometry B^ and Bg 
vanish both on the equator and on the axis, while B r van- 
ishes only on the axis and at the opening angle, where the 
largest field intensity exists. 

One remarkable prediction of these models is that the 
magnetic field at a given radius varies dramatically with 
angle (and hence probably with inclination to the line of 
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Fig. 5. Virial-isothermal low mass solution (a — 
—0.3,9 = 0.41): magnetic field maximum (in Gauss) as a 
function of the distance to the central object (in AU) for 
different angles (9 = 7,15,30,45,60). The hashed region 
corresponds to the maximum magnetic field for different 
self-similar index (a = -0.45, -0.3, -0.1,Af « 1 M Q ). 



100 Gauss at 1 AU (e.g. Hughes |1999| ), as shown in Fig.| 
There is also a weaker dependence on the self-similar index 
as is also shown in Figfi These predicted magnetic field 
strengths are surprisingly large. This is nevertheless con- 
sistent with existing (rare) observations and reinforce the 
idea that magnetic fields pla y a m ajor role around young 
stella r objects (Donati et al. 1997 , Guenther & Emerson 
19961) . 
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3.3. Massive Protostars 

To adapt the model for high mass protostars, we must 
find solutions with higher self-similar density (according 
to Eq. |l^) using the method outlined in Sect. [O] and 
adjust the other variables, mainly the magnetic field in 
fact, until the boundary conditions are satisfied. For ra- 
diative high mass protostars, we may speculate that we 
may not need a magnetic field as strong as in the low mass 
case to either launch or collimate the outflows since strong 
pressure gradients are prod uced b y radiative heati ng (see 
also Shepherd & Churchwell |i"996; , Churchwell ^99^) . How- 
ever the virial-isothermal model that we give below has a 
stronger magnetic field. 

3.3.1. The Virial-isothermal Case 

Fig.p shows the self-similar variables for a virial- 

isntriprmal massivp nhippt. Tlnis snlntinn ic alon pnrrmlpfplv 
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Fig. 6. Massive protostar virial-isothermal solution: Vari- 
ations with angle of the self-similar variables (related 
to velocity (u r , uo, u$), density (p, and fi a ), magnetic 
field (y p , j/0 ) and Mach number (M ap , M a <f>) with a = 
—0.1,0 = 0.1 (solid lines). The poloidal component S p of 
the Poynting flux is plotted in the lowest right panel. A 
low mass solution (M sw 1 Mq) is plotted with dashed 
lines over the high mass case (M w 30 Mq). 



which the radial velocity changes sign) is the same in both 
cases, we observe that massive protostars produce faster 
(scaled) radial velocities throughout the positive range. 
There is also a more negative ug throughout this range 
which leads to an enhanced 'collimation', corresponding 
to an increased axial density of material. In addition, the 
peak rotation at the 'turning point' and the infall velocity 
just before this peak are both larger for the massive star. 
The smaller Alfvenic Mach numbers also show that the 
magnetic field is more important in this flow. This seems 
consistent with the ability to confine the higher physical 
temperature. 

Our high mass virial-isothermal protostar yields larger 
scaled velocities than does the low mass object. This is 
mainly due to the larger luminosity in the high mass case. 
Sufficiently close to the axis both the low (1 Mq) and high 
(30 Mq) mass objects can produce 300 km a -1 . However 
as noted previously it is not clear that streamlines at such 
small angles are part of the present flow. 

3.3.2. The Radiative Case 

Massive molecular outflows present large bolometric lumi- 

nnRit.v and ttiprpfnrp radiation nrnnanlv rtlavs an irrmnr- 



Fig. 7. Massive protostar radiative solution: Variations 
with angle of the self-similar variables related to velocity 
(u r , ug, u^), density (p and /i a ), radiation field (f r , fg), 
Mach number (M ap ) and temperature (O) for a = —0.2. 
M w 5 Mq. 



ture are no longer constant with angle and the radiation 
flux is not purely radial. The index a.f in Eq. |7] is chosen 
to be negative (—0.1 as in FH1) and so simulates a slight 
radiation loss. The numerical method is similar to that 
used above. The self-similar radiation flux f(8) is fixed at 
dmin to be zero , and the non-zero value at 9 max measures 
the radiation 'loss' from the self-similar region. 

Fig.p] shows a solution for the radiative case where 
a = —0.2, a = and 6 = 2. The temperature O is at its 
maximum (0.5) in the jet region and decreases towards the 
equator. The components of the radiation flux are much 
smaller than in FH1. Thus the required heating is less im- 
portant than in FH1. 



Fig. 3.3.2 presents two illustrative radiative solutions 
for low (M (w 0.3 Mq) and high mass (M w 7 Mq) cases. 
Streamlines in the poloidal plane are overplotted on con- 
tours of the hydrogen number density in the upper part of 
the figure, while the intensity of the radiation field is rep- 
resented at the bottom. Fig. 3.3.2 clearly shows the saddle 
type singularity of the central point. We see that a Keple- 
rian disc naturally appears as part of the global solution 
as does the collimated outflow. 

[[Caption of the figure that is given in gif format be- 
cause of its size] Streamlines and density contours of the 
hydrogen number density (in logarithmic scale) for low 
and high mass radiative solutions with a — —0.2. The 
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The radiation field is smaller in the massive case in 



Fig. 3. 3. 2. This is remarkable since one should expect sig- 



nificant radiative heating surrounding massive protostars. 
However, we have chosen to have the same components 
of the self-similar velocity for both cases, as well as the 
same self-similar temperature and radiation flux as input 



parameters. Therefore Fig. 3.3.2 shows that the necessary 
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radiation field has to be larger in low mass outflows in or- 
der to get the same scaled velocity as in the massive case. 
Particularly in the latter case, the global shape of the ra- 
diation field i s rath er similar to the solution obtained by 
Madej et al. ( 1987 ), where an anisotropic radiation field 
was produced by scattering in thick accretion discs. It 
was shown in their article that collimation was clearly dis- 
played by the radiation field in the polar region. We also 
find, in the massive case, that the radiation field increases 
away from the equator. On the other hand the low mass 
case almost shows an almost spherical symmetry, except 
on the axis. 

Radiative heating increases the pressure gradient push- 
ing outwards and therefore peak and asymptotic velocities 
along streamlines are larger as represented in Fig.pL This 
figure shows variation of total velocity with integrated 
distance along a streamline for both radiative and virial- 
isothermal cases for the same input parameters. We begin 
the integration in the equatorial region at about 10 4 AU 
from the protostar and integrate out to 10 6 AU in the axial 
region, where we attain a final angle of 0.1 radian (see FH1 
for comparison). Velocity is maximum at the turning point 
(u r = 0), which is the closest point to the source, and is 
almost completely toroidal there. The inclusion of radia- 
tive heating provides a net acceleration of material which 
clearly increases the asymptotic axial velocity. Moreover 
the asymmetry in this graph indicates that the angular 
variation dominates the radial decrease as dictated by the 
self-similar forms. In fact this asymmetry is the clearest 
indication of the production of high velocity outflow from 
the relatively slow (but massive) infall. 

Global differences in the solutions between virial 
isothermal and radiative cases remain as described in FH1. 
The temperature is maximum in the jet region and de- 
creases towards the equator. The radial component of the 
radiation flux decreases monotically in the same direction 
and its ^-component increases, being always positive. 

In order to show the influence of the mass on solutions, 
we present variations of the total velocity with the inte- 
grated distance along streamlines as defined in the previ- 
ous plot for different masses in Fig|| For all the curves, the 
input parameters are kept the same (a = —0.2, = 0.5) 
except for the mass of the central object and therefore the 
self-similar density \i. Even the initial components of ve- 
locity and magnetic field remain the same. Only the start- 

inc distance nf intpcrratirtn fYnm tine r-pntral nhipr't varies si-v innnl naramptpr^ (Vplnn'tv and macmetif 1 fipld spalpd 



Fig. 8. Variations of total velocity with integrated dis- 
tance along a streamline for radiative (solid) and virial- 
isothermal (dashed) cases (a — —0.2, Q a xis = 0.5). The 
streamline are integrated from a radius of about 10 4 AU in 
the equatorial region out to 10 6 AU in the axial region, at a 
final angle of 0.1. The central object mass is 5 Mq. Arrows 
denote the direction taken by the flow along streamline. 



tiplied by approximately a factor 5 as material flows from 
the equatorial region to the axial region. The efficiency 
of heating as a driving mechanism grows with the input 
values of the temperature O and of the components of 
the self-similar radiation flux /. Indeed, a factor of 10 or 
more can be obtained just by increasing these parameters. 
Moreover, the figure shows that the highest velocity com- 
ponents (« 100 km s _1 ) are most tightly collimated along 
the axis (9 = 0.02). Thus the most massive protostars pro- 
duce the fastest flows with maximum radial velocities in 
the axial region. 

We conclude that radiative heating combined with 
Poynting flux driving provides an efficient mechanism for 
producing high velocity outflows when the opacity is dom- 
inated by dust. 

3.4. Parameter Study 

[[Caption of the figure that is given in gif format because 
of its size] Monte Carlo shooting: Plots of maxima of u^ 
at the turning point (u r = 0) with respect to the maxima 
of ue for different values of the temperature parameter 
(O=0.01;0.16;0.26;0.36). (a = -0.2)] 

We use a Monte Carlo exploration of our parameter 
space in order to study general trends in the model. The 
self-similar index and the temperature are fixed and the 
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opening angles are found at the lower left of the band for 
each temperature, that is with maximal poloidal collima- 
tion flow and the smallest possible maximum rotation. 

It is also found that larger opening angles are associ- 
ated with smaller magnetic fields. The corresponding re- 
sults are not reported here, but in fact the variation of the 
opening angle is rather large for a relatively small mag- 
netic field variation (a reduction of 25% in the magnetic 
field magnitude can widen the opening angle from 35° to 
60°). If the magnetic field varies secularly as the evolution 
proceeds, then a sequence of our models could be regarded 
as a series of 'snapshots' of the protostellar evolution. This 
evolutionary sequence would show that the opening angle 
increases with time, as the magnetic field becomes less im- 
portant. This would be consistent with the notion that the 
youngest outflows are generally the most collimated. In- 
deed, this has been observationally verified by Velusamy 



Fig. 9. Variations of total velocity with integrated dis- 
tance along a streamline for protostars of different masses. 
Solid lines represent a central object of 10 solar masses 
while dotted and dashed respectively stands for 5 and 1 
M Q . Two sets of solutions are represented for two stream- 
lines at different polar angles {6 — 0.02; 0.1). 

exists a broad range of solutions that possess favorable 
characteristics. Only the most significant results are re- 
ported here. 

We plot the maxima of u^ as functions of maxima 
of ug for various values of the self-similar temperature 
(6=0.01;0.16;0.26;0.36), and for a given index {a = -0.2). 
These maxima occur at the turning point where the ra- 
dial velocity vanishes. Each point shown here represents 
a 'good' solution. Rotation is directly measured by u^, 
and the tendency for material to go towards the axis, 
i.e. the collimation, is related to ug. We find, as a global 
trend, that when (ug) max decreases, {u^max remains al- 
most constant ((u^,) max « Q.2'b{ug) rnax + constant) un- 
til {ug) max reaches a threshold value. Then (u0) ma:E de- 
creases faster than (ug) ma x, which is what was anticipated 
analytically for super- Alfvenic flow when a = 1/4 (see ap- 
pendix B). In addition, there is a limiting (u<p) max for each 
(ug)max at a given temperature, below which solutions 
are not found. When increases, which physically cor- 
responds to larger pressure gradients, (u^) max decreases 
while (ug)max becomes less negative. If the magnitude of 
the latter quantity is taken as a measure of the collima- 
tion, we see that rotation and collimation decrease together 
with increasing temperature and therefore probably with 
the bolometric luminosity of the central object. 

It is found that solutions are less influenced by density 
(and hence mass) than by temperature. Nevertheless the 

HistriVintinn nf nnpninff ancrlp iq not fnnnH t.n phantrp its the two terms in Fn. 



and Langer (1998), who provide evidence that outflows 
widen in time. Moreover when the magnetic field reaches 
a lower threshold in our model, the solution becomes a 
pure outflow. Therefore, as obtained with our model, if 
the opening angle broadens sufficiently, it may ultimately 
cut off the accretion supplying new material for the out- 
flow. 

The angular dependence of the self-similar variables 
are represented in Fig.nft The only parameter that varies 
is the self-similar temperature. The figure (left panel) il- 
lustrates the fact that collimation decreases with temper- 
ature as already mentioned. Moreover variations of the 
self-similar temperature give rise to different density pro- 
files in the infalling 'disc' region, as seen on the right hand 
panel of Fig.[l0| The density decreases close to the equa- 
tor for the largest values of the temperature with a max- 
imum well away from the equator. Such solutions accrete 
so rapidly in the equatorial plane that they have a reduced 
density relative to their rotationally supported surround- 
ings. These solutions resemble the special case a = 1/4 
as discussed in appendix B. On the other hand the tem- 
perature has little effect on the density profile in the axial 
region, although it broadens the central 'funnel' substan- 
tially. 

We note that solutions for cold flows (O = 0) have 
also been obtained. In such a case, the circulation is pow- 
ered only by Poynting flux since our Bernoulli integral 
(Eq. |32|) applies. We do not expect asymmetric velocities 
at the two infinities on the stream lines. For the same set 
of parameters, the total velocity in the 'interacting' re- 
gion is smaller for cold flows and the fast rotating region 
(where u^ is maximum) is narrower relative to finite pres- 
sure flows. Moreover the Alfvenic Mach numbers almost 
reach unity at the peak rotation point, showing that the 
flow is only slightly super- Alfvenic there, and that in fact 
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Fig. 10. The angular dependence of ue (left panel) and of density fj, (right panel) for various self-similar temperature 
parameters (with a = —0.2). Solutions with 9 = 0; 0.3; 0.8 are represented with solid lines. Two consecutive solutions 
are separated by a change of 0.1 in temperature. 
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Fig. 11. Dependences on the self-similar index a — 
—0.1, —0.2, —0.3 (with = 0.3). Radial (upper panel) and 
toroidal (lower panel) components of the self-similar ve- 
locity are plotted as functions of the angle. 



are plotted for three values of the self-similar index but 
with the same self-similar temperature. For more negative 
self-similar index one gets dramatically smaller 0- and <f>- 
velocities (recall that a = —1/2 yields pure radial accre- 



4. Observational Consequences of the Model: 
Synthetic 13 CO Spectral Lines and Maps. 



In this Section, we compute u CO (J=l 



emis- 
sion lines for several of the outflow solutions discussed 



0) 



in Sect. 3.1. Specifically, we provide examples of spectra, 
channel maps, position-velocity diagrams, and intensity- 
velocity diagrams for one example of each of the following: 
the low mass solution of Sect. 3.2, the radiative high mass 
solution of Sect. [3.3.2 , and the virial- isothermal high mass 
solution of Sect. 3.3.1. The importance of this analysis is 



that it allows us to directly compare the observable fea- 
tures of our models with real observational data. However, 
we cannot fully explore the observational consequences of 
our models in the present paper. Our purpose in this Sec- 
tion will be mainly to outline the qualitative features of 
our models, leaving a complete exploration of the param- 
eter space for the second paper in our series. 

The method that we follow is essentially the same as 
that outlined in FH2, except that the code used there 
has been substantially improved in a number of respects. 
Firstly, we are now able to generate results with much 
higher spectral and spatial resolution. Secondly, although 
we compute spectra on a grid of pencil beams through the 
outflow source, we convolve our spectra with a Gaussian 
telescope beam to more accurately simulate observational 
results. FH2, on the other hand, did not perform this con- 
volution, which resulted in many spuriously sharp features 
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Fig. 12. Normalized spectra are shown for the low mass 
solution on a 5 x 5 grid of map positions that are evenly 
spaced across the mapped region, which is 20 x 40 r in 
size. The abscissa is in units of km s~ l . The source is 
located at the center. 



has no significant effect on the spectra or maps, except for 
in the lowest, and hence least interesting, velocity chan- 
nels. Since we are primarily interested in emission at rela- 
tively high velocities, corresponding to the outflow, emis- 
sion or absorption by relatively slow moving background 
gas is of no consequence. As in FH2, the primary limita- 
tion of our analysis is that we assume local thermodynamic 
equilibrium. Although this is probably not strictly true, we 
do expect the 13 CO level populations to be approximately 
in equilibrium with the local temperature provided that 
the optical depth is at least moderately high. Therefore, 
we hope to capture at least the character of the emission, 
if not the precise intensity. 

4-1- Low Mass Solutions 

We have assumed a mass of 1 M Q for the low mass solu- 
tion, which implies a fiducial radius of r — 1400 AU using 
Eq. O. We have also chosen a inclination angle of 45° be- 
tween the symmetry axis of the outflow and the plane of 
the sky. With the mass and r determined, the density, 
temperature, and velocity are determined at any position 
by the dimensionless parameters fi(9) and 0((9), and the 
radial scalings given by Eqs. U|j]. Therefore, by construct- 
ing a line of sight through the cloud, we can compute 
the emission (expressed as a radiation temperature) as a 
function of velocity parallel to the line of sight; examples 
of these spectra are shown in Fig.[12| We have computed 
spectra on a 61 x 120 grid of map positions across the 
outflow. 





-0.5 0.5 

Fig. 13. Channel maps for the low mass solution at an 
inclination angle of 45°. In units of km 8 , the velocity 
channels are as follows for panels a) through f): a) 5.5 < 
\v\ < 6.5, b) 4.5 < \v\ < 5.5, c) 3.5 < \v\ < 4.5, d) 
2.5 < \v\ < 3.5, e) 1.5 < |v| < 2.5, f) 0.5 < |v| < 1.5. The 
upper part of the outflow is the red-shifted side, while the 
lower part is blue-shifted. Map coordinates are in units of 
10 4 AU. 




Velocity Magnitude [km s 



Fig. 14. Intensity-velocity diagram for the low mass solu- 
tion. The dot-dashed line represents a power-law fit with 
an index of —3.5. 
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Fig. 15. Position- velocity diagram for the low mass solu- 
tion. 



relatively weak emission in the wings, as is the case for 
real outflows. Spectra may contain either a single wing or 
both red and blue-shifted wings depending on whether a 
particular line of sight crosses only a single lobe of the out- 
flow or both. There are also many cases where line wings 
are absent, but well-defined "shoulders" are present on the 
low velocity dominated spectra. 

In the equatorial regions of the maps, double-peaked 
spectra are found with relatively low velocity peaks. These 
spectra are mainly due to the radial infall that dominates 
in the equatorial regions. It is unlikely that the double- 
peaked profiles are due to rotation, since radial infall ve- 
locities exceed the rotational velocities at all angles 9 ex- 
cept for very near the radial turning point of the outflow. 
Furthermore, many of the spectra have a slight asymmetry 
in which the blue-shifted peak shows stronger emission. 
This is suggestive of the well-known outflow asymmetry, 
which is due to the more efficient self-absorption of blue- 
shifted emission by overlying material moving towards the 



observer (Shepherd & Churchwcll 1996, Gregersen et al 



The outflow is apparent at all velocities shown in our 
maps. We obtain a relatively wide outflow "cone" at the 
lowest velocities (0.5 to 1.5 km s" 1 ) and very well col- 
limated jet-like emission apparent at velocities between 
1.5 and 3.5 km s — 1 . At higher velocities, the jet-like fea- 
ture breaks up into outflow "spots" that move away from 
the central object as the velocity increases. The most im- 
portant feature of these channel maps is that the opening 
angle of the outflow gradually decreases as the magnitude 
of the velocity increases. As discussed in FH2, this effect is 
entirely due to the angular velocity sorting of the outflow 
solutions. Our models always have the property that the 
outflow velocity increases towards the axis of symmetry. 
Thus, our models naturally produce outflows in which most 
of the material is poorly collimated and moves at relatively 
low velocities, but the fastest jet-like components are very 
well collimated towards the axis of symmetry. This prop- 
erty has, in fact, been observed in a number of outflow 



sources (Bachiller et al. 199C, Guilloteau et al. 1992, Gueth 
et al |1996D . 

Several authors have noted that molecular outflows 
seem to be characterized by a power-law dependence of 
total integrated intensity as a function of velocity (Ro- 



driguez et al. 1982 Masson & Chernin 1992, Chandler et 
al. 1996 ). In Fig.|14|, we have summed all of the spectra in 
order to show how the intensity is distributed with veloc- 
ity. The solid line in the figure represents the blue-shifted 
emission, while the dashed line represents the red-shifted 
emission. We note that there are no important differences 
due to the somewhat low optical depth of our solution. 
We find that the intensity / is well fit by a power law 
I ol v 
0.2 km s 

found for several real outflows , and the index is in rea- 
sonable agreement with the available data (Cabrit et al 
1996, and Richer et al. 1998). It remains to be seen, how- 



35 over a vj velocities greater than approximately 
. This power-law behavior has, in fact, been 



ever, how sensitive the power-law index is to the parame- 
ters of our model. This important issue will be resolved in 
the second paper in this series, where we more completely 
explore the line emission of our models. 

that the emission 



An important feature of Fig. 14 



is 



1997, Zhou et al. 1996). The effect is slight in this case 
since the low mass solution, with a 1 M@ central ob- 
ject, produces a peak optical depth of only ss 0.2 along- 
most lines of sight. The massive solutions, discussed in 
Sect. 3.3.2 and 3.3.1 , produces much higher optical depths 
and leads to a much more pronounced asymmetry. 

Fig.O shows a set of channel maps for the low mass so- 
lution. To make these synthetic maps, we have integrated 
the emission over 1 km s _1 wide channels from —6.5 to 
6.5 km s _1 . We have deliberately excluded the emission 
with \v\ < 1 km s , since we find that the lowest velocity 

pmissirtn is psspntiallv fpatiirplpss Wp liavp pnrwnlvprl t.Vip tial pmissinn at. lnw vplnpit.ips alnncr t.Vip pnt.irp put Tlnp 



falls very rapidly with increasing velocity. For example, 
the intensity at even a rather modest velocity of 3 km s _1 
is only 0.01% of the peak intensity. At some velocity, the 
emission will fall below the noise threshold that is always 
present in real observations. Past this velocity, the emis- 
sion shown in Fig.O would be undetectable. 

In FigJIq, we show a velocity-position diagram for a 
cut along the outflow axis. We have actually computed 
spectra along several cuts parallel to the outflow axis and 
convolved them with a Gaussian beam to produce a single 
velocity-position diagram. We find that there is substan- 
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Fig. 16. Normalized spectra are shown for the radiative 
solution on a 5 x 5 grid of map positions that are evenly 
spaced across the mapped region, which is 10 x 20 r in 
size. The abscissa is in units of km s _1 . 

thus, wc find that the most intense emission occurs at low 
velocities at all map positions. 

The outflow is apparent in Fig.na by the emission ex- 
tending out to approximately ±8 km s _1 on either side 
of the outflow. Moreover, the general appearance of the 
position-velocity diagram is suggestive of a globally ac- 
celerated flow since higher velocities generally appear at 
greater distances from the central source. How can this 
be reconciled with the fact that velocity decreases with 
radius in our model? It was shown in FH2 that the appar- 
ent acceleration is due to the unique velocity sorting that 
is predicted by our model. The radial velocity always in- 
creases with decreasing polar angle near the outflow axis. 
Our cut crosses progressively smaller angles with increas- 
ing distance from the central object, and therefore encoun- 
ters higher velocity material. 

The ragged appearance of the highest velocity contours 
is entirely due to the limitations of our numerical proce- 
dures at high velocities. In our model, the highest velocity 
components are extremely localized near the central star 
and the outflow axis. Thus, whether or not a pencil beam 
passes through such a component is a matter of chance 
when the spatial extent falls below the grid spacing. This 
apparently happens when |v| > 5 km s . 

4-2. Massive Solutions 

4.2.1. The Radiative Case 

It is useful to compare the 13 CO emission of the radiative 



d) 



2 
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b) 



e) 



w 




-1 1 

Fig. 17. Channel maps for the radiative solution at an 
inclination angle of 45°. In units of km s , the velocity 
channels are as follows for panels a) through f): a) 11 < 
M < 13, b) 9 < |v| < 11, c) 7 < |v| < 9, d) 5 < |v| < 7, e) 
3 < |u| < 5, f) 1 < |v| < 3. The upper part of the outflow 
is the red-shifted side, while the lower part is blue-shifted. 
Map coordinates are in units of 10 4 AU. 



mass. Since the mass is a free parameter, we have tuned 
it to achieve realistic optical depths of order unity. The 
mass is larger than for the low mass solution discussed 



solution discussed in Sect. 3.3.2, with the low mass solu- 
tion Hisr'nQHPr] flhnvp We havp assiimprl a ma^ of 1 A/f^ 
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Fig. 19. Position-velocity diagram for the radiative solu- 
tion. 



Fig. 20. Normalized spectra are shown for the massive 
solution on a 5 x 5 grid of map positions that are evenly 
spaced across the mapped region, which is 2.5 x 5 r in 
size. The abscissa is in units of km s _1 . 



above, but the effects of radiative heating are more likely 
to become important in outflows surrounding more mas- 
sive protostars, in any case. The only other parameter that 
we have changed is the FWHM of the Gaussian beam that 
we convolve with our solution; here, we can use a slightly 
larger and more realistic 10 arcsec beam without smear- 
ing out the details of the solution. As in Sect. 4.1, we have 
computed spectra on a 61 x 120 grid of map positions 
assumed a 45° inclination angle. 

Qualitatively, we find a great deal of similarity with 
the low mass solution. The line profiles shown in Fig. 16 
have similar structure, and the outflow shown in Fig. 17 
is collimated to a similar degree. The main difference is 
the velocity of the flow. In the radiative case, we find that 
there is significant emission out to at least 10 km s _1 , and 
we find some low intensity emission out to 15 km s^ 1 ; this 
is most clearly illustrated by the position velocity diagram 
shown in Fig. |l9|. Examining the intensity- velocity diagram 
shown in Fig.|l8|, we find that the intensity / ex w -36 for 
velocities greater than about 0.5 km s _1 . It is striking 
that nearly the same power law index should be obtained 
for both the low mass and radiative solutions. Admittedly, 
we do not fully understand the reasons for this similarity 
in behavior at present. We also doubt that this sort of 
behavior is a generic feature of our models. Nevertheless, 
we are encouraged that at least some of our models are 
capable of reproducing such realistic observed properties. 

4.2.2. The Virial-Isothermal Case 

Wp have alsn ertTrmiiterl flip ^^dCl emission nf the hitrh 



fore, we have used a mass of 100 M Q and a corresponding 
fiducial radius of r = 44000 AU to compute the spectra. 
Even using such a large central mass, we find that the peak 
optical depth in the line centre is very high, with values 
exceeding 5 in many positions. Such high optical depths 
cause real difficulty for the present version of our code, 
since any given line of sight encounters the "photosphere" , 
where the optical depth at some velocity is approximately 
unity, rather abruptly. Thus, a great deal of care was taken 
to ensure that the solution was well-sampled along each 
line of sight. We note that the corresponding increase in 
computational time required us to decrease the spatial res- 
olution of our maps. Here, we employ a 41 x 80 grid, as 
opposed to the 61 x 120 grid used elsewhere. The inclina- 
tion angle is 45°, as in the previous two cases. 

The spectra shown in Fig.EtM are much more complex 
than the spectra of the low mass or radiative solutions. 
Many show multiple sharp peaks and very extended wing 
emission out to well past 20 km s _1 (as in Shepherd 
& Churchwell 1996, Gregersen et al. 1997 for CO maps 



of molecular outflows, and Doeleman et al. 1999 for SiO 



masers) . We are quite certain that the multiple peaks are 
real, and not numerical artifacts, since we have verified 
that their location and appearance is independent of the 
spectral resolution and the degree to which a solution is 
sampled along a line of sight. It remains unclear to us in 
detail why such jagged peaks should be present, although 
it must depend on the variation of optical depth, density 
and velocity field within the 'beam'. Thus behavior should 
be studied as a function of inclination to the line of sight. 

Desnit.p the nrmsnal armearpmee nf the srteetra the 
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Fig. 21. Channel maps for the massive solution at an 
inclination angle of 45°. In units of km s -1 , the veloc- 
ity channels are as follows for panels a) through f): a) 
30 < M < 35, b) 25 < M < 30, c) 20 < \v\ < 25, d) 
15 < \v\ < 20, c) 10 < \v\ < 15, f) 5 < M < 10. The 
upper part of the outflow is the red-shifted side, while the 
lower part is blue-shifted. Map coordinates are in units of 
10 4 AU. 

Fig]23 shows significant emission at most positions out to 
approximately 20 km s . We obtain emission at higher 
velocities (up to approximately 40 km s , but the emis- 
sion is restricted to the central map positions nearest the 
protostar. The intensity-velocity diagram (Fig.E2) again 
shows a power law behavior, but the slope is different than 
that obtained in the previous two cases; here we find that 
/ oc v , which is still inside the r ange a llowed by the 



Velocity Magnitude [km s ] 

Fig. 22. Intensity-velocity diagram for the massive solu- 
tion. The dot-dashed line represents a power-law fit with 
an index of —4.1. 
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avail able observations (Cabrit et al. |1996| , Richer et al 
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Fig. 23. 

tion. 



-10 10 

Velocity (km/s) 

Position-velocity diagram for the massive solu- 



5. Discussion 

5.1. The Alfvenic Point 

All the solutions presented previously were completely 
super-Alfvenic. Here we report a trans-Alfvenic solution 
in order to show why they are generally not preferred. 
In such a solution, the outflow size is reduced and the 
ultimate velocities are smaller. Therefore this type of so- 
lution is less efficient in producing high velocity outflows. 
Entirely sub- Alfvenic solutions can also be found but they 

admit nnln rrnrp nrrrptinn T^io' 941 sVinw« t.VlP Hpn«itv rtrn- 



critical Alfvenic density. In this case, the Alfvenic Mach 
number is larger than unity. The second case (B) has al- 
most the critical density at the turning point, and conse- 
quently the Alfvenic Mach number almost reaches unity 
at this point. However the solution always remains super- 
Alfvenic. Finally the lower panel (case C) displays a trans- 
Alfvenic solution showing that the size of the outflow has 
been reduced drastically. Therefore the most interesting 
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ing flux that could assist the outflow as in the present pa- 
per. The inclusion of this effect allows us to obtain faster 
and more collimated outflows. The path followed by ma- 
terial also differs from FH1 as shown in Fig.E3, where 
streamlines from FH1 and the present models are plotted. 
All the solutions start from approximately the same initial 
position in the equatorial region with the same set of pa- 
rameters. Streamlines from the present model get closer 
to the source than the FH1 solution, for both radiative 
and virial-isothermal cases. Thus, in the present model, 
the flow passes relatively close to the central mass. The 
limiting speed is about the escape speed from the central 
mass (and the jet speed) of a few hundred km s _1 . We 
note that the most collimated of the the three solutions is 
the radiative case. 



10 
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Fig. 24. We show the density /x and the Alfvenic density 
\i a for really super Alfvenic (upper panel), slightly super 
Alfvenic in the region of naught radial velocity (middle 
panel) and critical (lower panel) solutions. From top to 
bottom Poynting flux and radial velocity are getting less 
and less important (a = —0.2,9 = 0.41). 
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Fig. 25. Comparison of streamlines without Poynting 
flux (dashed line) and including Poynting flux (solid line) 
starting from the same position in the equatorial plane. 
(a = -0.3,6 = 0.41) 



5.2. Comparison with FH1 



5.3. Stability Analysis 

The steady configurations presented in this article should 
develop strong shocks in the axial region due to the strong 
collimation created by magnetic and pressure forces. 
Moreover current-carrying jets, as in the present model, 
are liable to Kelvin- Helmholtz (KH), Pressure Driven 
(PD) and magnetic instabilities driven by the electrical 
current. These so called current driven (CD) instabilities 
have been studied only recently in the context of astro- 
physical jets (Appl |1996| , Appl et al. |l999j ). We use the 
results of the latter linear stability analysis in order to 
get a rough estimate of CDI and KHI in our magnetic 
configuration. Growth rates T of the different modes show 
that the CD helical mode (m = 1) dominates the CD 
pinching mode (m = 0) and that CD and KH insta- 
bilities become comparable when the Mach number ap- 
proaches unity. This is the case in the present self-similar 
model. The CD instability (with a typical wavelength 
Acd ~ ^Routfiow, Rout flow being the characteristic size of 
the outflow) should be located on a current sheet where 
reconnection and particle acceleration might occur. KH 
pinching modes (Xkh « 3 x R ou tfiow) should dominate 
KH helical modes, and consequently give rise to internal 
shocks in the outflow. All these results require numerical 
simulations that we plan to study in another paper. 



6. Conclusions 

In this paper we have presented a model based on the r- 
self-similarity assumption applied to the basic equations of 
ideal axisymmetric and stationary MHD, including Poynt- 
ing flux. Detailed comparisons with the observations have 
been computed and characteristic scales of the problem 
have been given as functions of source properties. The lu- 
minosity needed for radiative heating is smaller even while 

Hrivina' faster outflows than in thp nrpviniK moHel fT^Hl 1 
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surround them, although the jets may be partly due to 
central activity not included in our model. 

The most massive protostars produce the fastest flows 
with maximum radial velocities in the axial region. Ra- 
diative heating produces faster outflows compared to the 
virial-isothermal case, for both low and high mass objects. 
Larger opening angles are associated with smaller mag- 
netic fields. Consequently, a gradual evolutionary loss of 
magnetic flux may result in outflows that widen as they 
age. 

Synthetic spectral lines from 13 CO (J=l — > 0) al- 
low direct comparison with observational results via chan- 
nel maps, maps of total emission, position-velocity and 
intensity-velocity diagrams. The model reproduces well 
observational features. Due to internal instabilities in the 
most collimated 'jet' part of the flow, the time evolution 
of the steady model should give rise to regularly spaced 
knots, and possibly excitation, ionization and non-thermal 
particle acceleration due to field annihilation and shock 
dissipation. Thus the model in its current state of devel- 
opment shows that radiative heating combined with Poynt- 
ing flux driving is efficient in producing high velocity out- 
flows when the characteristic luminosity of the forming 
star (as deduced from observations) is used. More exten- 
sive searches in parameter space, the detailed fit to spe- 
cific sources, and time dependent modeling all remain to 
be done. Although we do not discuss the regions excluded 
from the self-similar region of the flow, the current cor- 
respondence to observational features is sufficiently exact 
that we expect simultaneous infall/outflow ('circulation') 
to be an essential part of any realistic model. 
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Appendix A: The Equations 

The mass flux conservation equation remains the same as 
in FH1: 



(1 +2a)fXU r + -T—p;-j7, (tiugsm9) = 0, 
sin dti 



(A.l) 



while the self-similar variable y is simply replaced by its 
poloidal component in the following set of equations: 
Magnetic Flux conservation: 



(a + 5/4) u r 1 d fug sin 







y v sin 9 d9 \ y p 

Radial component of momentum equation: 



(A.2) 



i_fL 2 ^ e + i + ^^,o (A.3) 



conservation: 
1 du <t> , , 
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d$ 



= (A.4) 



Faraday's law plus zero comoving electric field: 



do 



d , 
UQUr + UQUe—m 

du 



1 
Vp 



-component of momentum equation: 



U r UQ 
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du r 



M„J d9 
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y P M ap 
d<d Odfi 
~dl + ~jid9 



de 



(A.6) 



To treat the radiative heating we either hold constant 
(virial-isothermal case) as discussed above, or we include 
the equations of radiative diffusion (radiative case) as in 
FH1. The equations (17), (18) and (22) of FH1 continue 
to apply provided that the optical depth given above is 
>2/3. 

Appendix B: The First integrals 



Equations ( A.l ) and (A2) together yield the first integral 
H(ue sin 6) W^yp^ 



<7i/4tt, 



(B.l) 



where y p replaces y in FH1. Moreover directly from 
Eq. A.l and the poloidal stream-line equation dr/rd6 — 
u r /ug we have the stream-line integral 



r {1+2a) iiu e sin ( 



m- 



(B.2) 



From this we see immediately that a > —1/2 if we are 
to have quadrupolar stream-lines for which r — > oo 
as ug — > 0. A second stream- line integral follows from 



Eq. A.2 and the stream-line equation, but it is not inde- 
pendent of Eqs. B.l and |B.2J since it follows by eliminating 



fi between these two. The limit a < 1/4 follows immedi- 
ately from Eq. B.l if we require finite /i, ug = and odd 
symmetry in the magnetic field at the equator. The inte- 
gral then implies implies that y p — > oo and so the field 
passes through zero in the equatorial plane. This conclu- 
sion is avoided on the polar axis since the radial velocity 
is free to become very large there. This choice has the ad- 
ditional merit that the equator is super Alfvenic as well 
as the axis. 
Using Eq. 



together with Faraday's law (AJj) we 
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Note that the constant q 2 is very interesting because it 
can be used to measure the strength of the electric field 
which is in fact 



E 



GM 



cr„ 




1 
Vp 



1 

!Jd> 



(«0 x u p ).(BA) 



The integral is not present in FH1 since the electric field 
is everywhere zero there. 

In general we can not find an angular momentum inte- 
gral explicitly (this is related to the difficulty of expressing 
scale invariance in an action principle; B. Gaffet, private 
communication) but in the special case of similarity index 
a = 1/4 such an explicit integral exists in the form (cf 
FH1): 



ug sin ( 
M 2 



q 3 (u<f, sin 9) 1 



1 



M ap M a< f, 



(B.5) 



Although this is a very special case not normally related 
to our numerical solutions (the mass density is constant 
in radius so that eventually the dominance of the central 
mass is broken as r increases) , it is instructive to consider 



it further. Eq. B.l becomes simply /i = {qi/Air)y, v l so that 



Ml v = qiy p and M^ = q x y\ly v . Then Eqs. |J and |[5 
can be solved together for y p and y^ to give 



1 _ 9293 sin 2 6 + qiqzu^ue sin 2 I 



Vp 



qzu^ug : 



U V< 



and 



1 _ Q.i<liU% sin 2 6 - q 2 

y<P 93"! sin 2 6* + u&Ue 



(B.6) 



(B.7) 



We therefore observe that near the equator y p can be in- 
finite only if Ucf, — > faster than ug . This feature is shared 
with our numerical solutions, but of course the equator 
is strictly outside the domain of the solution. From the 
second equation we see that y<f, tends to zero at the equa- 
tor under the same conditions. If one were to insist that 



u,j, 7^ at the equator, then Eq. B.5 shows that we need 
M a< pM ap = 1 there, or by the above Mach number defi- 
nitions y$ = l/<7i. Our solutions for y p and y,p show that 
this is only possible if q 2 — 0, that is with zero Poynting 
flux. Thus it seems from our analysis of this special case 
(however we always find u^ = at the equator in our 
solutions) that the presence of a non-zero Poynting flux 
is not compatible with a Keplerian equatorial disc. The 
material there is falling radially towards the star. This is 
a possible non-linear end state for an instability that cou- 
ples disc rotation to AlfVen waves propagating out of the 
plane (Shu et al. EL9941 Dendy et al. |1998|). 
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